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Since the breakthroughs in 2005 which have led to long term stable solutions of the binary black 
hole problem in numerical relativity, much progress has been made. 1 present here a short summary 
of the state of the field, including the capabilities of numerical relativity codes, recent physical 
results obtained from simulations, and improvements to the methods used to evolve and analyse 
binary black hole spacetimes. 
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I. INTRODUCTION 



With the current generation of gravitational wave in- 
terferometers (LIGO [il-[3, Virgo GEO (B-i] and 
TAMA [^) reaching their design sensitivity, and the 
next generation (Advanced LIGO [l^ and LISA [Ulil) 
planned to come online in the next few years, the accu- 
rate modelling of gravitational wave sources and the con- 
struction of template waveforms to be used in matched 
filtering of detector data has become urgent. One of the 
strongest expected sources of gravitational waves is the 
inspiral and merger of a binary black hole (BBH) system. 
Additionally, the modelling of BBH mergers has applica- 
tions in astrophysics, where it can be used to learn about 
possible growth methods for supcrmassive black holes, 
galaxy evolution and globular clusters. BBH mergers oc- 
cur in the regime of very strong gravity, and whilst the 
early inspiral can be described accurately with approxi- 
mate analytic post-Newtonian (PN) methods, to model 
the late inspiral and merger requires the numerical solu- 
tion of the Einstein equations using expensive and com- 
plex supercomputer simulations in numerical relativity 
(NR). Though the analytic framework and associated nu- 
merical simulations have been in development since the 
1960s, it was only in 2005 that simulations of multiple 
orbits of a BBH system became possible. Since those 
first simulations |13l4l5| , the quality of the computations 
has improved significantly, and it is now possible to simu- 
late enough orbits with a sufficiently high accuracy that 
the resulting waveforms can be compared with PN re- 
sults. NR waveforms can be attached to PN waveforms 
to construct /i?/&rirf waveforms, fitted with phenomenolog- 
ical models, or used to tune Effectivc-One-Body (EOB) 
waveforms. 

In this brief status report, I review recent developments 
in the field of numerical relativity. As such, this work is 
not intended to be a comprehensive discussion of the his- 
tory of the subject. For more details on earlier work, 
particularly associated with waveforms for gravitational 
wave (GW) data analysis, see Ref. [l^. Several works 
have appeared as this manuscript was being prepared, 
and this demonstrates the healthy level of activity in the 
field. This article is structured as follows. In Sec. [Hi I 
give a brief overview of the current capabilities of NR 
codes, including the longest waveforms and highest mass 



ratios studied so far. I also summarise the results of 
a validation study (the Samurai project) which verified 
the consistency of the waveforms produced by the several 
codes in the community, and conclude the section with 
a discussion of some work on improvements to computa- 
tional efficiency. 

In Sec. mil I discuss new physics results which have 
been obtained using NR simulations, including compar- 
isons with PN and EOB models. I then discuss some ap- 
plications of NR to gravitational wave data analysis, in- 
cluding work on injecting numerical waveforms into data 
analysis pipelines (the NINJA project), the construction 
of phcnomenological waveform templates, and studies of 
NR results applied to the detection of gravitational waves 
and the estimation of the parameters of their sources. I 
then discuss work that has been done on modelling the 
state of the final merged black hole as a function of the 
parameters of the initial black holes, as well as work on 
hyperbolic (non-circular) and high velocity BBH config- 
urations. 

In Sec. IIVI I present recent work on methods used in 
NR for BBH simulations. Proposed new techniques for 
generating improved initial data are reviewed, as well 
as work which has been done on boundary conditions. I 
then discuss studies involving the computation of gravita- 
tional waves from numerical simulations, including anal- 
yses of finite radius effects and extrapolation procedures, 
the asymptotic falloff of the Weyl scalars used for mea- 
suring spacetime content including gravitational radia- 
tion, and the first application of Cauchy-Characteristic 
Extraction (CCE) to produce unambiguous waveforms 
at future null infinity for BBH systems. I also discuss 
the development of new computational infrastructure for 
using non-Cartesian grids with finite differencing BBH 
codes, as well as improvements to the methods used for 
spectral BBH evolutions. 

In Sec|Vl I list the main challenges facing the NR com- 
munity today, and present some closing remarks. 



II. STATE OF THE ART 

There are now several groups around the world with 
codes capable of performing BBH simulations [H, [TBI . 
[T7l - [26j . There are two analytic forms of the Einstein 
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equations (BSSN [27l - [29| and Generalised Harmonic [ij, 
Si mi), and various numerical methods (high order finite 
differencing, pseudo-spectral, adaptive mesh refinement, 
multi-block) in use in the community. 

A. Consistency 

The gravitational waveforms computed by these codes 
for a given set of initial data should be independent of 
these differences as they represent a physical observable, 
so a comparison of the results of the different codes yields 
a strong test of consistency. The Samurai [s^l project was 
a collaborative effort to compare the waveforms from five 
different codes for the last 4 orbits and merger of a binary 
of equal-mass, non-spinning black holes in circular orbit. 
The codes were BAM CCATIE [H], Hahndol [ISiIll, 
MayaKranc ^ and SpEC [U. The first four use the 
BSSN formulation of the Einstein equations, though the 
coordinate conditions used are slightly different for each 
code. SpEC uses the Generalised Harmonic formulation. 
CCATIE and MayaKranc both use the same computational 
framework (Cactus/Carpet [s^ - ls?! ). but the evolution 
and analysis codes are separate. The other codes used 
in the Samurai comparison are independent. SpEC uses 
a pseudo-spectral evolution scheme, whereas the others 
use finite differencing methods. Despite these differences, 
it was found that the waveforms all agreed with each 
other within their published error estimates. In gravita- 
tional wave data analysis, the quantity usually used to 
make comparisons between waveforms is the match, Ai 
[ssf . This is a normalised frequency-domain scalar prod- 
uct of the waveforms, weighted by the detector noise, in 
the range zero to one, where a match of = 1 indi- 
cates two waveforms which are seen as the same by a 
given detector. Wc also refer to the mismatch 1 — M. 
of two waveforms. The mismatch between the Samu- 
rai waveforms was less than 10""^ for enhanced LIGO, 
Advanced LIGO, Virgo and Advanced Virgo, suggesting 
that the accuracy and consistency of the waveforms is suf- 
ficient for detection purposes. We should note here that, 
since only the last eight cycles of the waveform were con- 
sidered, this result applies only to those systems where 
these cycles arc the only ones in the detector bandwidth; 
i.e. those systems above a certain total mass. Also, only 
the dominant I = 2,m = 2 multipolar mode was consid- 
ered, approximately corresponding to a system observed 
from a direction perpendicular to the orbital plane. The 
impact of the inclusion of higher modes and longer wave- 
forms (e.g. constructed by joining NR to PN waveforms) 
in such a comparison has not yet been studied. 



B. Longest Waveform 

The longest NR BBH waveform so far produced lasts 
for 15 orbits and includes the merger and ringdown 
phases, and is described in Boyle et al. [s^ and Scheel 



et al. [40[. This waveform, from an equal mass binary of 
non-spinning black holes, was generated using the SpEC 
code and, due to its length and quoted accuracy, has 
been used in a number of studies comparing NR and PN 
results [4T| - |46| . This code uses a pseudo-spectral numer- 
ical method and multiple coordinate patches to cover the 
evolution volume and, in Szilagyi et al. j47| . new tech- 
niques were introduced which allow simulation of the full 
inspiral, merger and ringdown without the fine-tuning of 
parameters that was previously required. As a result of 
this work, simulations of BBH mergers with mass ratios 
q = mi/m2 > 1 (where mi > m2 are the masses of the 
individual holes) of g = 2 and dimensionless spins up 
to 0.4 are now possible with the SpEC code, and these 
are presented in Ref. [47| . Additionally, simulations with 
dimensionless spins of 0.44 anti-aligncd with the orbital 
angular momentum a re p resented in Chu et al. (48| . In 
a talk by H. Pfeiffer [4^, a series of long unequal mass 
simulations performed with the SpEC code was presented, 
with 15 orbits up to q = A and 8 orbits up to g = 6. 



C. Highest Mass Ratio 

For a 3D numerical relativity code, one of the most 
challenging problems is the simulation of black holes of 
very different masses. Unequal masses are expected to 
be commonplace astrophysically, hence waveforms from 
these systems are essential. One early reason for the 
interest in unequal mass BBH systems, besides the use 
of NR waveforms for GW science, was the fact that the 
mass asymmetry leads to an asymmetry in the linear mo- 
mentum emitted in the gravitational waves around the 
merger, leading the final black hole to recoil out of the 
initial zero momentum frame. This recoil (also known as 
a kick or rocket effect), especially when the black holes 
are spinning, can have important consequences for astro- 
physics [53 • The highest mass ratio that has been sim- 
ulated is q = 10. The first such simulation, reported in 
Gonzales et al. [5l| , consisted of three orbits of a binary 
of non-spinning black holes as well as the merger and 
ringdown of the final black hole to Kerr. The result was 
computationally very expensive to achieve, but it veri- 
fied the prediction of the empirical recoil formula js^ - is^ 
in a range which had previously not been tested. The 
reason for the increase in computational expense is that, 
for a fixed total mass M = mi + m2, the gravitational 
wavelength remains approximately constant with varying 
q, but the length and time scale required to resolve the 
smaller hole scales approximately with q. With adaptive 
mesh refinement as used in Ref. [5l|, for large q, the com- 
putational cost in CPU hours is q times higher than for 
an equal mass simulation. One technical problem which 
arises when simulating unequal mass BBH systems relat- 
ing to the coordinate conditions used in many codes has 
been studied recently in Miiller et al. [5^. Due to the 
computational expense of simulating high mass ratios, it 
can be desirable to propagate the gravitational waves in 
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a separate evolution. In Lousto et al. (56[, a BBH simula- 
tion with a mass ratio of 1:10 was presented, along with 
a computationally inexpensive perturbative evolution of 
the gravitational waves modelled using point particles as 
the sources, whose locations were determined from the 
coordinate tracks of the black holes in the numerical sim- 
ulation. 



D. Highest spin 

Mathematically, Kerr black holes have a maximum di- 
mensionless spin of 1, and there is a good probability 
that highly spinning black holes exist in nature 57 59j. 
Current NR techniques for simulating BBH systems are 
limited in the maximum spin which can be specified in the 
initial data, both for theoretical and numerical reasons. 
The common initial data formulations have mathematical 
limits on the spins which can be achieved (see Sec. IIV Al 
for a discussion of techniques for constructing initial data 
with high spins). Numerically, black holes with high spins 
require more resolution and hence higher computational 
resources to achieve a given accuracy. Typically, black 
holes with dimensionless spins as high as 0.6-0.8 can be 
evolved with only a moderate increase in computational 
cost over the non-spinning case. The cost of simulat- 
ing higher spins increases dramatically however. For the 
Bowen-York initial data used by the majority of BBH 
codes, the theoretical maximum has almost been reached 
in Dain et al. [gO], where ~ 7.5 orbits of black holes with 
dimensionless spins 0.92 are evolved. This is the highest 
spin for the inspiral-merger-ringdown of a BBH system 
so far simulated in NR. 



E. Lowest Eccentricity 

BBH systems visible to GW detectors are usually ex- 
pected to have circularised to an orbit of eccentricity zero 
due to the emission of gravitational radiation [6l| . How- 
ever, it is not known how to rigorously choose BBH ini- 
tial data parameters so that the resulting orbit is quasi- 
circular. The lowest eccentricities have been obtained 
using a method presented in PfeifFer et al. [g^]- In this 
method, the first few orbits of a BBH evolution are per- 
formed and the initial data parameters adjusted in an 
iterative procedure based on the measured eccentricity 
of the resulting evolution. With this method, eccentric- 
ities as low as e ~ 5 X 10~^ have been produced. This 
procedure must be repeated for every initial data set con- 
sidered, though the computational cost is mitigated by 
the fact that the evolutions do not need a very high res- 
olution. It should be noted that there is no fully general 
relativistic definition of eccentricity; definitions can be 
used which are based on the coordinate or proper sep- 
aration of the black hole horizons, and these definitions 
are used in Ref. [g^- It would also be possible, though 
more difficult, to use the properties of the gravitational 



radiation, such as its amplitude or frequency. 



F. Computational Performance 

As higher mass ratios and longer simulations arc re- 
quired, the performance of NR codes becomes more im- 
portant. The XiRel project [H, [13] was started in or- 
der to improve the performance of the publicly available 
Carpet [So . [37j adaptive mesh refinement infrastructure. 
As a result of recent work. Carpet now scales efficiently 
up to 2048 processing cores, and as it is used for BBH 
simulations at AEI, GaTech, LSU, RIT and UIUC, the 
improvements benefit simulations by all of these groups. 



III. NEW BINARY BLACK HOLE PHYSICS 
A. Comparisons with Analytic Models 

Due to the high computational cost of NR BBH sim- 
ulations, the early inspiral can instead be modelled us- 
ing approximate analytic techniques based on the post- 
Newtonian method (see Schafer 65| for a recent review), 
and the late inspiral and merger can be simulated with 
NR. PN results are accurate in the regime where the 
black holes are far apart and moving slowly. For GW de- 
tection purposes, the use of longer waveforms increases 
the range of masses of systems which can be detected, so 
including the inspiral is very important. Recent work has 
focused on combining the waveforms from analytic meth- 
ods with NR simulations to generate complete waveforms 
describing both the inspiral and merger. For a review of 
the first PN-NR comparisons see, for example, Hannam 

The Effective-One-Body (BOB) approach ^ (see 
Damour and Nagar [g^I for a recent review) is heavily 
based on the PN method, but uses additional techniques 
to model the merger phase, which is usually inaccessible 
to PN methods. There are several unknown parameters 
in the various FOB methods, and these need to be deter- 
mined by comparison with NR simulations. Two recent 
works, Damour et al. [4l| and Buonanno et al. [i^l, have 
used the SpEC 15 orbit equal mass non-spinning waveform 
[39I |40| to tune the parameters in the FOB formalism, 
adding to previous work on PN and FOB comparisons 
with this waveform in Refs. [i^, 13]. After this tuning, 
the FOB model agreed with the NR waveform within the 
quoted numerical error on the NR waveform. For binaries 
with total masses between 30 and 150 Mq (solar masses), 
there was a mismatch between NR and FOB of < 10^'^ 
for LIGO, Enhanced LIGO and Advanced LIGO. This 
work was very recently extended in Pan et al. [46| to in- 
clude the new SpEC waveform (48j incorporating spinning 
black holes. 

In Campanelli et al. [g^, a 9-orbit NR simulation of 
a configuration of spinning black holes of unequal mass 
with their spin axes misaligned with the orbital angular 
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momentum was presented. This is the most general BBH 
configuration that can be expected in nature, assuming 
that such systems will have evolved to have zero eccen- 
tricity by the time the NR portion of the waveform enters 
the frequency range of GW detectors. The mass ratio 
was q = 1.25 and the dimensionless spins of the black 
holes were 0.6 and 0.4. In this simulation, significant 
precession of the orbit out of the initial orbital plane was 
observed, in agreement with PN spin-orbit and spin-spin 
coupling predictions (see, e.g. Kidder for details). A 
comparison with 3.5 PN waveforms including spin effects 
was made, and the first 6 wave cycles (3 orbits) agreed 
within 1% (note that this is a time-domain inner-product 
measure, not the conventional frequency-domain match 
Ai usually used in gravitational wave data analysis) . 

The work of Peters [Hlj in the 1960s demonstrated that 
a binary system in an eccentric orbit which is inspiralling 
due to the emission of gravitational radiation will lose 
eccentricity and eventually circularise. Typical gravita- 
tional wave sources are expected to circularise before they 
can be observed, so studies of sources of gravitational ra- 
diation usually assume circular orbits. However, it has 
been suggested (Ml that there may be sources with 
non-negligible eccentricity which are detectable by Ad- 
vanced LIGO and LISA. In Hinder et al. [75|, the first 
comparison between an eccentric PN model and an NR 
simulation was performed. A 10 orbit simulation (20 GW 
cycles) with initial eccentricity e « 0.1 was presented, 
and an eccentric PN model was fitted to the resulting 
waveform. The two waveforms agreed within 0.1 radi- 
ans for the first 8 cycles, and by 5 cycles before merger, 
a dephasing of 0.8 radians had occurred. Two different 
PN models were compared, and the choice of PN expan- 
sion variable was found to have a significant effect on the 
accuracy of the PN approximation. 



B. Gravitational Wave Data Analysis 

Due to the weak nature of gravitational waves, and 
the difficulties in separating GW signals from local vibra- 
tions, the output of ground-based gravitational wave de- 
tectors is dominated by noise. As a result, sophisticated 
statistical algorithms must be used to extract physical 
signals corresponding to the detection of gravitational 
waves from BBH systems. These algorithms, employ- 
ing the technique of matched filtering, require accurate 
waveform templates corresponding to the sources to be 
detected. 

For low mass systems (those with a total mass of 
between 2 and 35 Mq), current GW detector data is 
searched for signals using PN template waveforms. This 
is because for high masses, the GW frequency of the 
merger is too high to be seen by the detectors, whereas 
for low masses, pure PN templates can be sufficient [76j . 
Between 25 and 100 Mq, waveforms generated using the 
EOB formalisrn (not including spin effects) and tuned to 
NR (EOBNR [11 [Z^-Izi) as weU as phcnomenological 



templates [80|, |8l| are currently used. For a recent sta- 
tus report on the searches of LIGO- Virgo data for signals 
from coalescing binaries, see Sengupta et al. [82|, and for 



more information on the use of NR waveforms in GW 
data analysis, see Hannam [l6j . 

Now that NR waveforms are available which include 
the merger phase, as well as a significant overlap with 
the regime of PN validity, it is important to test that the 
search algorithms work well with signals which include 
the merger. This work was started in the Numerical In- 
jection Analysis (NINJA) project [H-H^- NR waveforms 
from ten different groups were injected into the data anal- 
ysis pipeline software and the performance for detection 
and parameter estimation of these pipelines was charac- 
terised. The NR waveforms covered mass ratios up to 
q = 4, a range of spin configurations up to dimensionless 
spins of ~ 0.9, and up to 30 GW cycles. The waveforms 
were converted into time-series data that would be seen 
at the detector and Gaussian noise, designed to mimic the 
features of each detector, LIGO, Virgo and Geo600, was 
added to the signal. Various data analysis techniques 
were then used to attempt to detect the signals in the 
data, and in some cases to estimate the source param- 
eters. Overall, it was found that many of the current 
data analysis pipelines were able to detect the merger 
waveforms at the expected sensitivities, but that signif- 
icant work is needed to improve parameter estimation. 
Additionally, the NINJA project has provided a working 
model for collaboration and communication between the 
NR and DA communities; essential if the recent advances 
in NR are to contribute to GW science. The NINJA 2 
project |85| . which is currently underway, will incorpo- 



rate hybrid waveforms (i.e. PN and NR waveforms com- 
bined), and real detector noise. This will allow stronger 
statements to be made concerning the performance of the 
detector pipelines for real signals. 

In Reisswig et al. [s^, a study was made of the sig- 
nal to noise ratio (SNR) of a family of waveforms for 
BBH systems with spins aligned and anti-aligned with 
the orbital angular momentum. In this family, there is 
no precession of the orbital plane. It was found that 
the SNR increases with the projection of the total spin 
in the direction of the orbital angular momentum, and 
that if the spins are aligned with the orbital angular mo- 
mentum and maximal (dimensionless spin of 1), the SNR 
is three times as high as in the anti-aligned case, lead- 
ing to an increase in the event rate for current and ad- 
vanced LIGO and Virgo detectors of a factor of ^ 30. 
Hence it is likely that there will be a bias in detected 
binaries towards high spins aligned with the orbital an- 
gular momentum, as these have the largest SNRs. It was 
also determined that the match between different cases is 
controlled by the total spin; binaries with zero total spin 
(including the case with both black holes non-spinning) 
are indistinguishable from each other (indicated previ- 
ously in Vaishnav et al. [l^) for LIGO, whereas binaries 
with a nonzero total spin have identifiably different grav- 
itational wave signatures. This means that aligned spin 
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detection templates only need a single spin parameter, 
which reduces the dimensionality of the search parame- 
ter space. 

In Ajith et al. [STj . it is shown that if only non-spinning 
templates are used to search for signals from spinning 
systems, the event rate is reduced by up to 50%. To alle- 
viate this problem, it is necessary to include spin effects 
in the templates used in GW data searches. The authors 
constructed time-domain hybrid waveforms by matching 
waveforms from NR simulations covering at least 8 wave- 
form cycles to a PN approximant called TaylorTl [ssj . 
Mass ratios up to g = 3 and dimensionless spins up to 
0.85 aligned and anti-aligned with the orbital angular 
momentum were considered. Phenomenological Fourier- 
domain waveform models with free parameters were then 
constructed, and the parameters were fitted to the PN- 
NR hybrids. This resulted in a bank of template wave- 
forms with freely adjustable mass ratio and spin param- 
eters which can be used in a GW data-analysis pipeline, 
and should lead to higher detection rates for spinning 
BBH systems. 

Once a detection has been made, the next step is to 
determine the parameters of the source, for example the 
sky location, mass ratio, total mass, etc. For the space- 
based detector LISA, the low noise levels mean that de- 
tection of BBH signals is generally not a problem; the 
challenge is to accurately determine the source parame- 
ters. In Mc Williams et al. [s^, a comparison is made be- 
tween the accuracy of parameter estimation when using 
purely PN waveform models and when using combined 
PN and NR models. It is found that the inclusion of the 
NR merger portion for supermassive BBH systems (for 
mass ratios q < 10) reduces the uncertainty in the sky- 
location of such a binary by a factor of 3 to within 

10 arcminutes. 



C. Modelling the Final State 



the models. In Refs. [96|,|97[, the authors proposed for- 
mulae for properties of the final black hole (mass, spin, 
linear momentum) based on symmetry arguments and an 
expansion in the initial spin parameter. The final spin 
has more recently been add resse d using an alternative ap- 
proach (also see Refs. f98l - ll00l | for predicti ons b ased on 
point-particles). InBarausse and Rezzolla |l0l| . follow- 
ing work in Refs. jl02l4l05 |. a set of approximations and 
simplifying assumptions is combined with a large body 
of already-published NR data to produce a general for- 
mula for the spin magnitude and direction of the final 
black hole after a merger. Instead of basing the model 
on PN expressions around the time of the merger, the 
final spin is expressed in terms of the spins early in the 
inspiral, and it is argued that this is more relevant when 
making astrophysical predictions, since an astrophysicist 
is not interested in the detailed dynamics of the merger, 
but instead in the final state obtained from a set of ini- 
tial parameters when the black holes are far-separated. 
In Lousto et al. Il06l . PN evolutions and the empirical 
model of Ref. [9^| are combined to make a link in a sta- 
tistical sense between the early inspiral and pre-merger 
spins. 

While the above works have concentrated on deter- 
mining the final state of t he b lack hole using empirical 
models, in Lovelace et al. 107| . the dynamical momen- 



tum flow in the region near the black holes is studied. 
The method used is based on field theory on an aux- 
iliary spacetime in which linear momentum is defined, 
and f ollow s on from previous work which treated the PN 
case |l08j | . In this way, the linear momentum of the in- 
dividual black holes can be computed, albeit with some 
gauge dependence. It is found that, despite this gauge 
dependence, the linear momentum measured for a head- 
on collision of two black holes using this method agrees 
well between evolutions using different gauges, suggest- 
ing that this linear momentum measurement may have a 
physical interpretation. 



Knowledge of the mass, spin and linear momentum of 
the final black hole in a BBH merger has many applica- 
tions. One example is in the modelling of the evolution 
of cosmological black hole spins [s^. Another is in ex- 
plaining the presence or absence of massive black holes 
in the centre of galaxies, as the merger of galaxies and 
their associated central massive black holes may result in 
recoils which eject the final black hole [50l |. 

When it was predicted and then confirmed [9ll - l93| that 
very large recoils could be generated in BBH mergers 
with specific spin orientations, a large amount of effort 
went into exploring this effect due to the potential ap- 
plication to astrophysics. Recently, in Lousto et al. j94| . 
the mass, spin, and recoil velocity of the final merged 
black hole were modelled as functions of the spins just 
before the merger using expressions based on PN pre- 
dictions, but with PN coefficients replaced by those ob- 
tained from fits to NR data. This is an extension of 
previous work [ssl loil . l95j including additional terms in 



D. Beyond Circular Inspirals 

Due to its astrophysical importance, the BBH problem 
in NR is usually studied for quasi-circular configurations; 
i.e. those consisting of circular motion combined with a 
low inward radial velocity. However, there have recently 
been several works which study non-quasi-circular con- 
figurations. 

Hyperbolic BBH encounters can be thought of as or- 
bits of eccentricity e > 1. In a Newtonian system, such 
a configuration would result in scattering of one black 
hole off the potential of the other, but in full GR, for a 
sufficiently small impact parameter, the black holes be- 
come gravitationally bound due to the emission of en ergy 
through gravitational waves and merge quickly jl09| . In 
Healy et al. fTToj . it was found that hyperbolic encoun- 
ters of spinning black holes resulted in a significant linear 
momentum imparted to the final black hole due to the 
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asymmetry of the emitted radiation at the moment of the 
merger. Recoil velocities as high as 10000 km/s could be 
obtained, in contrast with the quasi-circular case where 
the highest recoil veloci ty ob tained is of the order of 3300 
km/s. In Healy et al. it was found that the final 

dimensionless spin could be as high as 0.98 when extrap- 
olated to the case of maximally spinning black holes. 

In particle physics, high energy collisions arc used to 
probe the fundamental properties of nature. As these 
collisions arc performed at higher and higher energies, it 
is expected that they will result in the formation of an 
event horizon, and that these scattering processes will 
be well-modelled by the high velocity scattering of black 
holes. In Sperhake et al. jll2| . the authors studied high 
velocity head-on black hole collisions. The simulations 
included initial velocities up to w = 0.94c and, when ex- 
trapolated to black holes moving at the speed of light, 
the energy radiated in the head-on collisi on w as deter- 
mined to be about 14%. In Shibata et al. [ll3l |. the au- 
thors studied non-head-on high velocity collisions. By 
varying the angle of approach of the black holes, it was 
determined that the impact parameter should be approx- 
imately b < 2.5GM/c^/v for a black hole to form, where 
M is the total black hole mass. This model was fitted 
from the simulated velocities from 0.6c to 0.9c . A similar 
study was performed in Sperhake et al. |114| where the 



result was confirmed for these velocities, but for higher 
velocity {v ~ 0.94), the formula in Ref. [ll3l | may have 
been an overestimate. 

Zoom-whirl orbits are characterised by quasi-circular 
motion (whirl), alternating with highly eccentric motion 
(zoom). For BBH systems, this phenomenon has been 
seen in extreme mass ratio inspir als a s well as PN mod- 
els, and in Pretorius and Khurana Il09l it was seen for the 
first time in NR. In Sperhake et al. [114l | , zoom- whirl be- 
haviour is confirmed for t he h igh velocity collisions stud- 
ied, and in Healy et al. jll5l |. the authors study zoom- 
whirl orbits for a variety of mass ratios and spins and 
show that they can occur without the need to fine-tune 
the initial data parameters. Furthermore, in Gold and 
Briigmann [ll6j |. low momentum zoom- whirl orbits, ex- 
pected to more closely reflect astrophysical scenarios, are 
simulated and the dependence of the energy radiated as 
a function of the initial angle of the black hole momenta 
is computed. They find that as the initial angle between 
the momenta and the coordinate line joining the two BHs 
increases from zero, there is a peak in the emitted energy 
at 47 degrees, followed by several local maxima. 



E. Electromagnetic Counterparts 

Several recent works have studied possible electromag- 
netic (EM) signatures of BBH mergers. If such signa- 
tures are visible to current or future instruments, this 
opens up the possibility of coincident detection via EM 
and GW observations of BBH inspiral events. Not only 
might this increase the dctectability of these events, but 



it might significantly increase the accuracy with which 
the parameters of the binary, such as the sky location 
and redshift, can be determined. 

In Palenzucla et al. [ml . the fully coupled 

Einstein-Maxwell system is evolved for an equal mass, 
non-spinning BBH system in the presence of an external 
magnetic field but otherwise in vacuum. The system is 
designed to model the effect of the magnetic field sourced 
by an accretion disk around a supermassive black hole bi- 
nary, where the matter from the accretion disk has been 
cleared from the immedi ate n eighbourhood of the black 
holes by binary torques [119| . It is found that the EM 
fields evolve in a way which follows the dynamics of the 
binary and so they can be used as tracers of the motion. 
Additionally, the EM energy flux is found to oscillate 
with a quarter the orbital period of the binary and the 
energy in the EM field is found to increase as the binary 
approaches merger. For the astrophysically realistic field 
strengths used here, the addition of the EM field has no 
effect on the binary d ynam ics or waveform. In foUowup 
work by Mosta et al. |12G| |. therefore, the EM fields are 
treated as test fields and their evolution does not feed 
back into the evolution of the geometry. Here, larger ini- 
tial black hole separations are used along with black hole 
spins aligned and anti-aligned with the orbital angular 
momentum. It is found that for the lowest multipolar 
modes, the EM radiation follows very closely the ampli- 
tude and phase evolution of the GW radiation, so that 
detecting either one gives direct information about the 
other. Furthermore, the EM energy fiux scales quadrati- 
cally with the black hole spin, and is 13 orders of magni- 
tude lower than the GW energy flux. Most importantly, 
the frequencies of the EM radiation are well outside the 
range of existing radio observations and hence direct ob- 
servation is highly unlikely. 

In an alternative picture of the late inspiral of a su- 
permassive black hole system, if the surrounding gaseous 
environment is sufficiently hot, as for example in the nu- 
clear regions of some low luminosity active galactic nuclei 
(AGNs), there could be gas present ne ar th e black holes 
all the w ay up to and through merger 12l| . In van Me- 
ter et al. |122| , first steps are taken towards adding a gas 
into their BBH simulations by computing geodesies of 
the spacetime, and considering these to represent fiows 
of test particles on the background geometry of the BBH 
solution. Differences in the collision and outfiow speeds 
of the test particles are observed between single black 
hole and BBH systems, and between spinning and non- 
spinning black holes in binarie s. In Bode et al. |l23j and, 
independently, in Farris et al. |l24| , a supermassive BBH 
system is simulated in the presence o f a g as modelled 
using full GR hydrodynamics. In Ref. |124| . the rate of 
accretion of the gas onto the black hole is computed, and 
it is found that the luminosity of the resulting EM radi- 
ation is enhanced significantly over that expected from 
a single black hole, and esti mate s of dctectability of this 
radiation are given. In Ref. [l23l | , correlations are found 
between the EM and GW emission. For the most massive 
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supermassive black hole systems that would be detectable 
by LISA, the EM luminosities are here found to be high 
enough to be accessible to observations out to a rcdshift 
of z w 1 using X-ray measurements. 



IV. NEW BINARY BLACK HOLE 
METHODOLOGY 



A. Initial Data 

Specification of initial data for a BBH evolution is a 
nontrivial task as this data must satisfy the Einstein con- 
straint equations. These equations are nonlinear so it is 
not in general possible to construct initial data by su- 
perposing two single black hole solutions. There are two 
main methods for constructing BBH initial d ata in use 
for NR evolutions; Bowen-York puncture data |l25l . Il26j | 
and a method based on the Conformal Thin Sandwich 
(CTS) decomposition of the equations involvin g excision 
of the black hole interiors from the domain [62], 127l - ll3Cll | . 
In both these cases, very roughly speaking, one chooses 
the desired characteristics of the black holes and then 
solves a set of elliptic equations numerically. Both of 
these methods are subject to several problems: incorrect 
initial radiation content, difficulty in controlling eccen- 
tricity, and difficulty in including very high spins. 

Radiation content: In the majority of current simula- 
tions, the initial data contains spurious radiation (also 
known as junk radiation) in the neighbourhood of the 
black holes, due to simplifications made when construct- 
ing the initial data (in particular the conformal-fiatness 
assumption). It also does not contain the radiation which 
would have resulted from the evolution of the binary from 
early times. The spurious radiation is visible in BBH 
waveforms only at early times, and is followed by the 
astrophysically expected waveform. Typically the initial 
part of the waveform must be discarded, and the spurious 
radiation can also cause unphysical numerical reflections 
from mesh refinement or outer boundaries. In Johnson- 
McDaniel et al. 



131| , the authors propose a new method 



of constructing initial data based on the PN approxima- 
tion. This initial data set is not conformally flat, and 
contains the radiation expected from the previous inspi- 
ral. One drawback of this method is that the Einstein 
constraint equations are s atisfied only to a given PN or- 
der. In Kelly et al. (l32j . another typ e of initial data, 
also constructed using the PN method '133^, is evolved. 
This initial data is also constructed to satisfy the con- 
straint equations to a given PN order, as with previous 
approaches 134j . The resulting evolution gives the ex- 
pected waveform from the start and contains significantly 
less spurious radiation than the evolution of comparable 
Bowen-York initial data. On the other hand, there is sig- 
nificantly higher orbital eccentricity (e ~ 0.1), as well as 
a variation in apparent horizon mass, perhaps caused by 
the constraint-violating nature of the initial data. 



Eccentricity: Recently, in Walther et al. [135j, a PN 
method was used for constructing low eccentricity initial 
data parameters, incorporating an EOB or Taylor Hamil- 
tonian with Pade or Taylor expanded flux. Thi s me thod 
was compared with previous PN methods l20l. 13611 and 
eccentricities comparable to those in Rcf. [l36j are ob- 
tained for equal-mass non-spinning systems. Addition- 
ally, low eccentricities arc obtained for systems with un- 
equal masses and spins for a small selection of spin con- 
figurations. 

High spins: Both the Bowen-York and CTS initial data 
approaches have intrinsic limits on the spins of the black 
holes which evolve from the initial data. This issue has 
been addressed recently in Lovelace et al. [l37j . where a 
comprehensive overview of the problems of high spin ini- 
tial data is also presented. Whilst the dimensionless spins 
of the black holes on the initial data slice for Bowen-York 
and CTS initial data arc found to be 0.9837 and 0.99 re- 
spectively, this spin quickly relaxes to about 0.93 after a 
short period of time (sec also Dain et al. [g^l where this 
was previously found for the Bowen-York case). A third 
type of initial data is considered, based on superposed 
Kerr-Schild solutions (with constraint solution), and this 
is shown to allow relaxed spins of above 0.99. Note that a 
modification to the standard puncture approach for head- 
on collisions can be u sed to produce high-spin puncture 
data as well |l38l Il39| , and for the single black hole case 
jl40l | , though this has not yet been studied for the orbit- 
ing case. 

Constraint solution: When setting up a configuration 
of black holes for numerical simulation with Bowen-York 
initial data, it is necessary to solve the Einstein con- 
straint equations, which are elliptic equations on the 
initial data slice. This can be computationally expen- 
sive, and the most commonly used solver code for this, 
uses a bi-spherical coordinate sys- 
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TwoPunctures 

tem which is specifically adapted to the case of two black 
holes. To increase fiexibility, for example to investigate 
the possi bility of simulating more than two black holes 
(see Ref . jl42l | for an alternative approach), in Bode et 
al. [mI an approximate analytic method was used to 
solve the constraint equations, based on Ref. |l44l |. It 
was found that the main effect of this approximation 
was to alter the eccentricity of the resulting evolutions, 
and when only the merger and ringdown are in-band, 
the mismatch between the waveforms with the exact and 
approximate constraint solution schemes was < 10~^. 

Modifications to puncture initial data: Several studies 
of the properties of commonly use d ini tial data have re- 
cently been performed. In Brown [l45l |. the propagation 
of modes in the puncture geometry is studied with spher- 
ically symmetric BSSN evolutions, and it is found that 
a potentially problematic numerical reflection of pertur- 
bations in the geometry does not occur due to a decou- 
pling of the modes of the system. When initial data is 
constructed in the puncture technique, the coordinate 
conditions used in the evolution cause the geometry to 
change during the very first portion of the evolution and 
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settle to th e so- called trumpet solution. In both Han - 
nam et al. jl46l | and Immerman and Baumgarte |l47j . 
initial data is constructed which is already in the trum- 
pet form . However, no significant difference is found in 
Ref. jl46l | between the evolution of trumpet and standard 
puncture data. In current BBH simulations, the evolu- 
tion is usually performed on spatial slices which reach 
to spatial infinity. However, gravitational radiation is 
properly defined only at future null infinity (see the dis- 
cussion in Sec. lIV C|) . An alternative is to construct slices 
which asymptote to null infinity rather than spatial infin- 
ity, and one way to do this is to use hyperboloidal slices. 
In Buchman et al. jl48l l , Bowen-York initial data includ- 
ing unequal masses, spins and boosts is generalised to 
such hyperboloidal slices. The evolution of this initial 
data should allow extraction of gravitational radiation at 
future null infinity, removing a common source of error 
in BBH simulations. 



B. Boundary Conditions 

Typical NR BBH simulations begin with a three- 
dimensional initial data set and evolve this forwards in 
time. Due to the difficulties of simulating an infinite 
spatial domain with finite numerical resources, a finite 
domain is usually simulated in combination with an arti- 
ficial boundary condition on the exterior of the domain. 
Ideally, these boundary conditions should result in a so- 
lution which is indistinguishable from an evolution with 
an infinite spatial domain. This can only be achieved ap- 
proximately, and the boundary conditions should have 
certain properties to give a useful solution. Firstly, 
they should not contaminate the solution with unphysi- 
cal gravitational radiation, either due to reflections of the 
waves generated by the binary, or due to radiation gen- 
erated by the boundary condition itself. Secondly, they 
should be constraint preserving to yield a result which is 
a solution to the Einstein equations. Thirdly, they should 
result in a well-posed initial boundary value problem. This 
is a mathematical property which is a necessary condi- 
tion for the numerical schemes used to be formally stable. 
These three properties are often only approximately sat- 
isfied, but the effects can be included in estimated error 
bars on the physical results of the simulation. 

Current evolutions using the SpEC code in the Gener- 
alised Harmonic formulation such as those in Refs. (39l . 
l40j use an outer boundary condition designed to sat- 
isfy the Einstein constraint equations as well as_to min- 
imise the incoming gravitational radiation [3l|, Il49l ll5C| , 
though these are not designed to be mathematically well- 
posed (a necessary condition for formal stability of the 
probl em u nder small perturbations). Recently, in Wini- 
cour [l5ll | , boundary conditions for the Generalised Har- 
monic formulation are presented which are constraint- 
preserving and of the Sommerfeld type, which limits re- 
flections of waves off the boundary. Additionally, the 
initial boundary value problem is well-posed. 



BSSN simulations generally apply an outgoing radia- 
tive boundary condition by modelling the solution at th e 
boundary as a radially propagating outgoing wave jl52l | . 
This approach is designed to reduce reflections of the 
outgoing radiation into the numerical grid which would 
contaminate the solution. However, it is neither math- 
ematically well-posed nor does it respect the Einstein 
constraint equations. It is expected that the errors in- 
troduced are smaller than the errors due to finite differ- 
encing, but this has not been comprehensively studied 
for long simulations. One approach is to place the outer 
boundary sufficiently far from the binary that the wave- 
forms measured at finite radius are out of causal con- 
tact with the outer boundary, and this was the approach 
taken in the early short simulations as well as in modern 
simulations for which this is sufficiently inexp ensiv e (for 
example those in [l53j |). In Niuicz and Sarbach [l54j . new 
boundary conditions for the BSSN evolution system are 
proposed which exactly satisfy the constraint equations, 
specify approximately zero incoming radiation, and are 
stable (the initial boundary value problem is well-posed) 
to linear order. It is hoped that these boundary con- 
ditions can be used in evolutions in future to avoid the 
above-mentioned problems. 



C. Gravitational Wave Extraction 

Mathematically, gravitational radiation from a BBH 
inspiral and merger is defined at future null infinity; i.e. 
the part of the spacctimc towards which all null rays, such 
as light or GWs, propagate. Current NR BBH codes are 
based on a Cauchy, 3-1-1, formalism, which means that 
initial data is specified on a (portion of a) 3-dimensional 
spacelike slice and then evolved forwards in time. This 
leads to a sequence of spacelike slices on which the solu- 
tion is known. GWs must be read off at a finite spatial 
radius, and this introduces an error in the waveform. An 
obvious solution is to compute the radiation at several 
radii and extrapolate the waveform to infinite radius as 
a function of retarded time from the source. It was shown 
in Hannam et al. Ref. |155| that this is not a trivial pro- 
cedure, and a detailed analysis and successful extrapo- 
lation of the inspiral portion was presented in Boyle et 
al. Ref. [3^ . A failure to extrapolate waveforms from the 
SpEC code can lead to a ph ase error of 0.5 radians [s^, or 
a mismatch of 10~^ |l56l |. if using waveforms extracted 
at r = 50 M only. 

A detailed study of two different techniques for pcr- 
formin g th is extrapolation is reported in Boyle and 
Mroue |l56l | . The first method involves extrapolating the 
amplitude and phase of the complex waveform at each in- 
stant of time using various orders of polynomial extrap- 
olation. Essentially, this is the technique which is now 
common in the literature, though in Ref. [l56l | various 
refinements are proposed which are shown to improve the 
quality of the ex trap olation. The second method, based 
on ideas in Ref. [l55j, involves expressing the amplitude 
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and time as a function of the phase and then performing 
the extrapolation. Error estimates are made, and the two 
methods are shown to agree within these error estimates. 
However, there were problems with the convergence of 
both methods after the merger which were not properly 
understood, and future work will investigate these. This 
study applies only to the waveforms from the SpEC code, 
and the same conclusions might not hold for waveforms 
produced by other codes. Generally, numerical relativists 
have found that it is necessary to be extremely careful 
with the extrapolation procedure, and, as pointed out in 
Ref. |l56| . there are still issues to be resolved. 

In PoUney et al. |l57l |. wave extraction techniques at 
finite radius are studied and it is shown that the conven- 
tional extrapolation of waveforms extracted at low radii 
is consistent with extracting at r = lOOOM. The use of 
such a large extraction radius is possible because of the 
constant angular resolution in the NR evolution code, as 
demonstrated for lower extraction radius in Refs. j39l[6^ . 
The outer boundary can be placed out of causal contact 
with the wave extraction sphere without incurring inor- 
dinate computational cost. The peeling property of the 
Weyl scalars is also studied, and while V's a-nd "02 de- 
cay with distance from the source with the expected ex- 
ponent, the situation for ipQ and ipi is more complicated, 
and different exponents are obtained. One suggested ex- 
planation for the differing exponents is coordinate and 
frame effects in the wave extraction procedure, but it 
is also pointed out that the peeling property might not 
be expected to hold even for single spinning Bowcn-York 
black hole spacetimes, and by assumption it would not 
then hold for binaries either. 

Instead of performing a traditional Cauchy evolution, 
where the solution is evolved along timclikc directions, it 
is possible to evolve along null directions. This is called 
characteristic evolution. This cannot be done straightfor- 
wardly in the neighbourhood of the black holes as it leads 
to caustics in the solution, but it can be done in regions 
sufficiently far from them. The Cauchy-Characteristic 
Extraction (CCE) method combines the two approaches. 
Metric data is read at a sufficiently large radius from the 
solution generated by a Cauchy code and is used as in- 
put data for a characteristic evolution. In this way, the 
waveform can be evolved all the way to future null infin- 
ity, removing ambiguities due to coordinate conditions 
or insufficient distance from th e source. This has been 
implemented in Reisswig et al. [l58l Il59l | and represents 
the first unambiguous determination of GWs at future 
null infinity for a NR BBH simulation. Note that the 
outer boundary condition of the Cauchy evolution can 
still contaminate the waveform, but in this work the outer 
boundaries were causally disconnected from the CCE ex- 
traction worldtube, avoiding this problem. It was found 
that the waveform agreed with the extrapolated finite- 
radius waves within 0.2% in amplitude and within 0.01 
radians in phase. It was also determined that the dif- 
ference between extrapolated and CCE waveforms was 
not significant for gravitational wave detection, but that 



for Advanced LIGO, Advanced VIRGO and LISA, these 
differences might be important for parameter estimation. 



D. Numerical Methods 

For numerical relativity evolution purposes, the four- 
dimensional BBH spacetime is usually foliated with 
three-dimensional spacclike slices. Each of these slices 
can be split into two regions with very different compu- 
tational requirements. The region around the black holes 
has a complicated geometry reflecting the shapes of the 
horizons, but the region far from the black holes con- 
sists only of gravitational radiation. Since the radiation 
propagates essentially radially, it requires a constant an- 
gular resolution to resolve it. The majority of BBH NR 
codes today use Cartesian-type coordinates everywhere 
in the grid. These have the advantage of simplicity, as 
only a single coordinate patch is required to cover the en- 
tire simulation domain. However, they are not efficient, 
as they lead to an increasing angular resolution with ra- 
dius. If box-in-box mesh refinement is employed, this 
inefficiency can be partially addressed, but only at the 
expense of signific antly reduced radial resolution. 

In Pazos et al. [l60l |. multiple locally Cartesian coor- 
dinate patches adapted to the black holes and exterior 
spacetime were used with finite differencing methods and 
Generalised Harmonic evolution equations to compute 
an equal mass non-spinning BBH inspiral using excision 
techniques. The scalability and convergence were good, 
though the simulation deve lope d instabilities before the 



merger. In PoUncy et al. |153l |. a similar patch struc 



ture was used for the exterior spacetime, but adaptive 
mesh refinement (AMR) with Cartesian grids was used 
for the region surrounding the black holes. The entire 
inspiral, merger and ringdown was simulated. In both 
these works, the fact that the exterior spacetime is re- 
solved using constant angular resolution means that the 
computational cost in this region scales linearly with the 
radius. This means that the outer boundary of the do- 
main can be placed out of causal contact with the region 
of physical interest for only moderate computational cost. 
This bypasses the problem that the correct outer bound- 
ary conditions for the BSSN formulation for the BBH 
problem are not known. 

While most NR codes use finite differencing methods 
with global Cartesian coordinates, the SpEC code uses 
multiple coordinate patches as well as spectral methods. 
For simple equations, these methods can be shown to 
be significantly more accurate (exponential rather than 
polynomial convergence) for a given computational cost, 
and indeed, the quoted accuracy and efficiency on the 
waveforms from the SpEC code is impressive. In Szilagyi 
et al. H^l , a new gauge condition is presented along with a 
new grid structure (touching grids with penalty boundary 
conditions, rather than the previously used overlapping 
grids) , as well as a new method of dynamically adapting 
the grids to the shapes of the black holes. With these new 
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techniques, SpEC is able to handle a variety of configura- 
tions without fine-tuning the simulation parameters on a 
case-by-case basis during the merger phase as was previ- 
ously required. Simulations of mass ratios of (? = 2 and 
dimensionless spins of 0.4 are presented. 



V. CHALLENGES AND CONCLUSIONS 

In this work, I have given a brief overview of the cur- 
rent state of BBH simulations in NR. It is now possible to 
simulate BBH systems over periods as long as 15 orbits 
for low mass ratios and moderate spins, or to simulate a 
smaller number of orbits for mass ratios up to 1:10 or di- 
mensionless spins up to 0.92. Pushing these simulations 
to higher numbers of orbits, mass ratios and spins can be 
achieved but at significant computational cost. For arbi- 
trary spin configurations, it is not known how many black 
hole orbits are needed from NR before a sufficiently accu- 
rate complete waveform using input from analytic meth- 
ods can be constructed. In order to span the parameter 
space of initial spins and mass ratios, many simulations 
will be needed. These are still very expensive to perform, 
taking weeks or months on up to a thousand CPUs for 
high mass ratios and large numbers of orbits. 

Approximately speaking, the computational cost of a 
BBH simulation scales with the mass ratio 9, based on the 
resolution requirement of the smaller black hole alone. 
For an accurate equal mass, non-spinning 8 orbit simu- 
lation which takes on the order of 200 000 CPU hours, 
a simple estimate shows that scaling this to a mass ratio 
of g = 10 would take 2 million CPU hours, and this is 
probably an underestimate, as the fact that the waveform 
gets longer with increasing q has not been considered, nor 
the fact that the higher resolution required around the 
smaller hole will not always come with perfect computa- 
tional scaling. Since many simulations for different spin 
configurations would be required to construct or tune 
template banks used in gravitational wave data analysis 
for detection and parameter estimation, this is clearly in- 
feasible (computer resources for NR BBH groups for an 
entire year are usually of the order of a few million CPU 
hours). 

There are also theoretical issues wh ich n eed to be bet- 
ter understood. As explored in Ref. |l56l |. the extrapo- 
lation of finite radius waveforms is not yet perfect, es- 
pecially in the neighbourhood of the merger. Cauchy- 
Characteristic Extraction for BBH systems |l58l . Il59j | 
should in principle solve this problem, but its use is not 
yet widespread. It is also not yet clear exactly how to 
construct accurate waveforms from PN and NR results 
when the PN waveforms come with their own intrinsic 
errors. It has been shown [t^ that most of the different 
PN approximants can be used equally well for detection 
with ground-based detectors of systems below a certain 
total mass (12 Mq), indicating that the errors in each of 
these approximations are not significant for the low fre- 
quencies for which such systems arc in-band. However, 



for higher masses, and for parameter estimation, the PN 
errors become significant, and assuming that they are 
zero when constructing PN-NR waveforms is not neces- 
sarily the best choice. 

The generation of the gravitational recoil around the 
BBH merger is not understood in precise terms. Sim- 
ple analytic models have been tested [H, |9ll, [U, UHl and 
shown to fit well with numerical data. However, the mod- 
els are a posteriori empirical fits to forms qualitatively in- 
spired by PN predictions in a region where the agreement 
with PN is poor. Also, the crucial identification of the 
orbital plane is gauge dependent. To get a robust, coordi- 
nate independent understanding of this recoil will require 
detailed understanding of the physics in the regi on near 
the black holes. The work of Lovelace et al. |l07l |. where 
the linear momentum in this region is co mput ed, may 
be a starting point for this. Previous work |l6l| has also 
analysed the features of the recoil, showing that a Newto- 
nian model containing some input from the simulations 
combined with a Kerr quasinormal mode expansion is 
able to reproduce the num erical results. There has also 
been recent work |162| . |163| combining PN wit h pe rturba- 
tion theory in the close-limit approximation jl64l | where 
the recoil for non-spinning black holes is computed and 
agrees moderately well with NR results. The extension of 
this work to black holes with spins should lead to greater 
insight into the generation of the recoil. 

There have been recent theoretical advances which are 
not yet routinely incorporated in numerical simulations. 
More astrophysic ally r ealis tic initial data prescriptions 
have be en g iven |l3ll . Il33| . and initial evolutions per- 
formed [l32| . It is not clear to what extent the con- 
straint violations present in these initial data sets will 
impact their usefulness, and this needs to be studied fur- 
ther. A method for constructing initial data for almost 
max imally spinning black holes has also been proposed 
[l37j | which goes beyond the capabilities of the Bowcn- 
York initial data currently used by many groups. Two 
possibilities for improving the physical accuracy of grav- 
itational waveforms from BBH simulations have rece ntly 
been explored: hyperboloidal slicing conditions |l48l |. in 
which future null infinity is approached in the compu- 
tational d omain, and Cauchy Characteristic Extraction 
158lll59| . where the waveforms are evolved to future null 
infinity using a separate null code. Whilst it is shown 
that the CCE waveforms are consistent with currently- 
used extrapolation techniques, this picture may change 
as simulations become more accurate. 

In summary, there are many interesting problems to be 
addressed in vacuum NR. In terms of physics, these in- 
clude gaining improved understanding of the strong field 
geometry around the black holes, and how this affects 
the final black hole produced in the merger. Computa- 
tionally, increased efficiency for long inspirals, very high 
spins, and higher mass ratios are the main challenges. Fi- 
nally, a strong interaction with the community of grav- 
itational wave data analysts will be necessary to fully 
leverage the exciting progress that has been made in NR 
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to date. 
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